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1.      INTRODUCTION 

The  purpose   of  the   FORTRAN  program  DFASUB  is   to  integrate    a 
system  of  equations — described  "by   a  mixture   of  ordinary  differential 
equations,  nonlinear  equations,    and  linear  equations — using  a  discrete 
variable  method.      The  equations   are  written  in  the   form 

f(y_»  y_'  ,   t)  +  Pv  =  0  (1.1) 

where  the  vectors  y_  and  y_'   =   dy_/dt   are   of  length  p,   v  is   of  length 
q-p ,  P  is    a  q  x  (q-p)    matrix,   t  is  the   independent  variable,    and  f  is 
a  vector   function   of   length  q. 

DFASUB  uses   the  variable  values   to   integrate   the   system  using 
a  multistep  method  designed  for  use  with  stiff  ordinary   differential 
equations  which   also  works  well  with  non-stiff  equations    and  which  will 
handle  nonlinear  equations.      This  paper  deals  with  both  the  theory  behind 
the  integration  method  and  the  program  itself.      Section  2  presents 
the  theory. 

At  present,   the  routine   calls   a  number  of  other  subroutines  which 
are   compiled  by  a  general -purpose   system   [l].      The   function   of  these  will 
be   described  so   they   can  be   replaced  by  user-supplied  FORTRAN   subroutines. 
The   form  of  the   internal  variables   of  the  program  is      described  in  Section   3 
to   allow  the  user  to  write   subroutines   equivalent  to   those   above    as    described 
in  that  section.      Section   h  details    calling  parameters  to   allow   a  user  to 
write  a   complete   integration  package  built   around  DFASUB.      Section   5   describes 
the   detailed  program  logic  of  DFASUB. 

This   introduction,   together  with  Sections    3  and  k,   serve   as    a 
user's    guide  to   DFASUB.      Reading  Section  2   is   sufficient  to  understand  the 
theory  behind  the  routine.      Someone  interested  in  the  programming  aspects 
of  the  routine   may  read  Sections  1,    3-5   to  understand  how  DFASUB  works. 


2.   THEORETICAL  BACKGROUND 

Suppose  that  the  values  of  y_  are  known  (approximately)  at  a 

number  of  equally  spaced  points  t   .  ,  i  =  l,2,...k,  where  t.jn  >  t.. 

n-i  l+l  1 

The  backward  differentiation   formula  gives   the   relation 

hK  =   ~   3^  ( Vn  +   ai^n-l  +-  ■  '+   Vn-k} 

where  h   =   t      ,    -  t.    >   0.       (The    coefficients    a.    and   B„    can  be    found  in 
i+l  l  l  0 

Gear    [2],    p.    217.)      If  this    is    substituted  in    (l.l)    for  t   =   t     we    get 


(2.1) 


F   ( 


)    =   t(y9    - 


n  ^-n'   -i. 


+   E,t)+Pv     =0 
n       n  n— n 


(2.2) 


where 


E      =    -  -=—  (a  y  +, 

n  h3Q        Pti-1 


+   \  ^n-k} 


is   known.      Hence,    (2.2)    is    a  system  of  q  equations    in   the   q   unknowns 
y      and  v   .      In   general,   they   are  nonlinear.      If  these   are   solved  by   i 
Newton-like   iteration,  we   get 


where 


^n,(m+l) 


v 


-n,(mKL) 


^,(m) 


-n,(m) 


r-1 


J    "  Fn(^n,(m)'  ^n,(m) 


3F 


J  = 


3(i,v) 


3f 


a°  £  i  p 


3y_       h6n    3v_'    , 


(2.3) 


(2.U) 


and  y      ,    »    and  v      ,    v    are   the    q  iterates    for   the    solution  y      and  v      of   (2.2) 
ti,[i)  -n,(m)  "hi  -n 

If   accurate   initial  guesses  y      ,-N    and  v      ,<s    are  used,   very   few  iterations 

-n,(0)  -n,(0) 

of   (2.3)    are  necessary.      In   fact,   accuracy  is  not  needed  in  v 


as  we 


-n,(0) 

show  below  that  the  iterates   y      ,_,  N    and  v      ,_,  s    are   independent  of  v      /_.v 

*4i,(l)  -n,(l)  -n,(0) 

if  round-off  errors    are   ignored.      To  see   this,    compare  the   first   iterate 


y      , N    and  v      ,_x   with  the   first   iterates 
•41,(0)  -n,(0) 


f      /    v    and  v      /    x    starting  with 

y      /_,  x    and  v      ,_  *    starting  with  y      ,  _  N    and  v      ,  n*  .      From  (2.3) 
«4i,(l)  -n,(l)  s  ^.(O)  -n,(0) 


z  = 


4,(1)    "Zn.CL) 


^,(1)   "  ^n,(i: 


o 


•4i9(0) 


_:4i,(0) 


-  J"1  P    (v 


-    V 


nv    n,(0)    '"  -4i,(0) 


) 


Hence , 


Jz  =  J 


^1,(0) 


^1,(0! 


-Pn(V(0)   V(0)} 


3f 


a0    31 


3y_      hSQ  3y_' 


P     -  P 
n  n 


o 

j4i,(0)    "  ^n,(0)_ 


=   0 


Since  we   assume  J   is   non -singular   (or  the  Newton   iteration  will  not  work)  , 

we  see  that   z_  =  0;   hence,  f      ,    v    =  y      /_  \    and  v      ,    <.    =  v      ,    s    are   independent 

of  v      /«\.      Note  that   the   only   condition   on  J   for  this   to  be  true   is  that 
-n,(0) 

the   right-hand  p   columns    of  J  be  exactly  P    .      The  left-hand  q-p   columns    can 
contain   anything,   so  we  need  only   approximate  J   in  those  positions. 


Initial   Guess    for  y      ,    >. 

Since  past   values   are  known,    a  good  initial   guess    can  be   found 
using  polynomial  extrapolation.      We   assume  that  we   know  y      .  ,      i    <_  i    <_  k. 
d  y"      .      (The  latter  was    found  at   the  last  step  or   could  be   evaluated 


an 


from  (2 .  $  with  n-1   replacing  n.)      With  those  values  we    can  use   a  Hermite 
interpolation   formula  in  the   form 


Since  equal  intervals    are  used,   the  a.    and  3      are   independent   of  n   and  h. 

The  Basic  Algorithm 

The  initial  guess   is    calculated  from  (2.5) — it   is    called  the 

predicted  value,      v      ,„N    is   set  to  v     _,  .      Equation    (2.3)    is   iterated 
*  -n,(0)  -n-1 

until  v      /    \    and  v      /    \    appear  to  have    converged.      (if  they   do  not,   the 
■Hi,(m)  -n,(m) 

step  size  h  is    reduced  as    discussed  later   and  the  process   is    repeated.) 

The  error  introduced  by  one  step  of   (2.1)    is  known    to  be  of 

the   form 

,k+l      ( k+1 ) 

k+1  K^J 

If  the  error  in  the  numerical   solution   is  a  sufficiently  smooth   function 

of  h   and  t,   h  y  (^-i)   +  '-'(h  can  ^e    computed  using  a  difference 

formula  where   £     is   not  necessarily  the  same  point  as    £ ,  but  both   £      and 

£    are   in  the  interval   (t     ,  ,   t    ).      This   enables   the  error  to  be   estimated 

n-k       n 

under  the   assumption  that  y  changes   slowly   over  the  interval. 

If  the  error  estimate   is    too  large,    the  result  is    rejected,    and 
a  smaller  step   is   tried.      Otherwise,   the  step  is    accepted  and  a  suitable 
step  size   and  order   (k)    is    chosen   for  the  next  step. 

Internal  Representation   of  the  Algorithm 

We   have   used  a  Nordsieck    [3]  vector  to   represent  past   values   of 
the   data  because   of  the  ease   of  order   changing  and  error  estimation.      We 

will   describe   this   vector  below  for  a  single  variable  y      and  its    associated 

B —  n 

past  values  y    ,...y  ,  .   This  allows  us  to  use  the  vector  notation 
n-1     n-k 


^=    [yn>  hyn'  K-l'-"yn-k+l] 

for  the  set   of  saved  information   concerning  y   at  t    .      (T  is   the  transpose 

operator.)      When  we  step  from  t  to  t    ,  we   use    (2.5)    to  predict  y      ,  _.N 

n-1  n  n,(0j 

from  y  At  the   same  time,  we   save  the  values   of  y     ,  . . .  .y     ,     ,    from 

J4i-1  n-1  n-k+1  . 

y        .     We   can   represent  this  process  by 


n,(0) 


n-1 


n-2 


n-k+1 


61     a2 


0        0 
0        0 


0        0 


0        0        0 


Vi 

\ 

0 

0 

0 

0 

0 

0 

i 

0 

.. 

*n-l 


(2.6) 


Next  we   iterate    (2.3).      To   do   this  we   must    substitute  y      ,    N    and  v      /    >, 

",n,(in)     -n,(m) 

into  F  (y  ,  v  ).   Thus,  from  (2.2)  we  must  evaluate 
n  n   n 

a 

F  (y  ,    s  ,  v  ,  v )  =  f(y  (    v,  -  — —  y   ,  x  +  E  ,  t  )  +  P  v  .    .      (2.7) 
n  n,(m)    n,(m)      Jn,(m)    h3   n,(m)    n   n     n  n,(m) 

We  will  normalize    (2.l)    so  that  a     =   -1.      Let  us  write 


0 

Jn,(m)  8Q     n,(m)  n 


=hyn,(m-l)   +   3^(yn,(m)    "  ^.(m-lp    if  m  >   ° 


=  +  37  [aiyn-l  +   Vn-2  +"-+  Vn-k  +  hBl  hyn-l] 


(2.8) 


~-t  [aiyn-l+--'+   Vn-k]      if  m=  ° 


Hence, 


hyn,(0)    =" 


ak"ak  Pl 

n-1  3Q     yn-k        BQ     yn-l 


=  yv         + . . .  +  y,  y     ,    +  6  ^  hy '    n 
'rn-1  ■      Vn-k  1  °n-l 

Thus,    (2.7)    can  be  written   as 

Fn(jrn,(m)>  vn,(m)>  =  f(yn,(m)'  hK,U)'  V  +  Vn,(») 

If  ve  vrite Aj(m)  =    [yn>(m).  ^j(m).  Vl yn-k+l]T'  ve  see  that 

we    can   combine    (2.9)   with    (2.6)    to   get 


(2.9) 


(2.10) 


^,(0) 


TT        31      a2 


Yl      61     Y2 


10        0 
0        0        0 

0        0        0 


ak-l  ak 


Yk-1     Yk 


0  0 

0  0 


Zn-1  =  %-l 


(2.11) 


This   is   called  the  predicted  value   of  y.      Now  we  note   from  (2.8)    that 

2n',(mM)    =  *n,(m)   +  ^(yn,(irH-l)    "  yn,(m)}  (2'12) 

1  T 

where    c  =    [l,  — ,   0,.  .  .  ,   0]      and  y      /    . -,  \    -  y      /    \    is   given  by   a  component 
—  ^n  n,(m+lj        'n,(mj 

of   (2.3). 

The   final   step  is  to  perform  the  transformation  to   a  Nordsieck 

vector.      The   approximations  y    ,   y'  ,  y      _,,... y     ,     .,    determine    a  unique  k-th 
^  Jn'   Jn'  °n-l'        Jn-k+l 

k    (k) 

degree  polynomial.      Let   its    scaled  derivatives   at  t     he  y    ,  hy'  ,...   h  y     '   /k! 

n  n  n  n 

There   is   a  unique  non-singular  k+1  by  k+1  matrix  Q  such  that  if 


^[V^>,yM2>/2 *<*w 


then 


z     = 
— n 


%L- 


With  this  we  restate    (2.1l)    as 


an 


d   (2.12)    as 


where 


and 


Vo 


^1,(0) 


z 
— n 


=  QBQ        z     ,    =  Az      _ 
-n-1         -n-1 


,(m+l)   "  QZn,(m+l) 

=  ^n,(m)   +  I  (^n,(m+l) 

I  =  Qc 


^n.dn) 


A  =   QBQ 


-1 
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0      1 

(2.13) 


(2.1U) 


The  fact   that   A  is    a  Pascal  triangle   matrix  follows    from  the  fact  that  the 
predictor  performs   a  k-th  order  exact  approximation    z      ,    v    to   z    .      The 

advantage   of  this    form  is   that   the  step   size    can  "be    changed  easily  by 

2  k 

multiplying  the  vector   z     by   C(a)   =   diag[l,   a,   a    ,...,   a    ]  where 


a  =  h        /h   ,  ..      The  order   can  be   reduced  by   dropping  a  column   from  A,   a 
new     old  J  **      ° 


row  from  A  and  z    ,    and  changing  l_.      The   order   can  be   increased  by   adding 


— n 


k+l    (k+l) 
a  row   and   column   to  A,    changing  &_,    and  estimating  h        y  /(k+l)!.      This 

is   estimated  by  the  last  element   of   (z      -   z      _.)/(k+l).       (This   is   also 
J  -n       — n-1 

used  in   computing  the   error  estimate.)      The   values   of  l_  can  be   found  in 
Gear    [2],   p.    217- 

The   choice  of  step  size    and  order  depends   on   several  things. 
First,   if  the  error  estimate   after  the   corrector  converges   is   too  large, 
the  step  size   is   reduced  to  h/k   and  the   step  is   taken   over;    if  this    fails 
to  work  after  three  such   attempts,   the   order  is   reduced  to   1  and  a  final 
attempt   is  made  before   giving  up. 

If  the    convergent   corrector  value    is  within   error  bounds  ,  then 
every  k  steps   at   order  k  the   step  size   at  the   original  order,   order  k-1, 
and  k+l  is    calculated  with  a  bias   toward  the  lowest   order   (with  the   least 
work  required)    to   get  the   largest  step  size.      This   is   only  done   every  k 
steps   due  to   stability   considerations    [h]. 

The  error  estimator 


h  "4i  M_n     ,i   ,,n-.i  ,i      n 


j_0     J      n_J  J      n~J 

is   the  local  truncation  error  for  the  method  in  use   if  a     =   -1.      If 
y(t)   has    a  continuous    (k+l)-th   derivative,   then 

Lh(y(t))  =  Ck+1y(k+1)U)  hk+1 

if  the  method  is   of  order  k   and  t-h   <   5   <  t.      We   can   define 


v*k+1) 

c. 


k+1       hk+1(k+l)! 

since   for  a  k-th  order  method  applied  to  the   (k+l)-th   order  polynomial 
t        ,  the  error  term  L   (t        )    is  known  and 

(tw)(ktl)  -  (k*l)l 
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3.      REPRESENTATION   OF  VARIABLES;    SUBROUTINES 

DFASUB  requires   several  variable   arrays    for  integrating  a 
system;   these   arrays    are   described  below.      The   remaining  arguments   in 
the    calling  sequence   are  briefly   described  following  those.      For 
notational    convenience   N   is    the  number  of  equations   q,   NY   is    the  number 
of  nonlinear  equations   p,   and  NL  =  N-NY.      The  notation  A(b)C  means 
"from  A  to   C  in   steps    of  B" . 

Y(J,I)    for  I   =   l(l)NY   and  J  =   l(l)T   contains   the    current  value 
of  the  differential  and  nonlinear  variables   in  its   first   row.      Each 
subsequent   row  J   contains   the    (J-l)-th   derivative   of  the  variable  y. 
multiplied  by  H**( J-l)/ ( J-l) ! .      Thus,  Y   contains    all  the   elements 
necessary  for  a  Taylor  series   expansion  without   requiring  the   derivatives 
to  be   multiplied  by  H   /k!    before   use.      Thus,   the  predictor  step    consists 
of  multiplying  Y(l,J)   by  the  Pascal  triangle  matrix  as    described  in 
Section  2. 

YL(l)    for  I   =   l(l)NL  holds  the  value  of  the   linear  variables, 
corresponding  to  v  in  Section  2. 

DY(l)   holds   the  result  of   calling  the  evaluation   subroutine 
DIFFUN(T,G,DY,Y,YL,HINV) ;    for  I   =  l(l)NY   it    contains   the   difference 
y!    -   f(y.)    for  y.    =  Y(l,l),  y!    =  Y(2,l)*HINV  where  HINV  is   the  inverse 
of  the  step   size  H.      For  I   =  NY+l(l)N  the   error  in  the   linear  equations 
stored  in  YL(l)    is   in  DY(l).      Thus,   DY(l)    is   essentially  the    correction 
to  the  values   in  Y(l,l)   and  YL(j)   before  the  Newton  iteration. 

SAVE(J,I)    for  J  =   1(1)7,  I   =   l(l)NY  holds   the   initial  values 
of  the  nonlinear  variables    at   the  beginning  of  an   integration  step.      If 
the  step   fails,   the   independent  variable   T(l)    is   restored  to   its   original 
values    and  Y(j,l)    takes    on  the  values   in  SAVE(J,l).      YLSV(l)    for 
I   =  l(l)NL   serves  the   same  purpose    for  the  linear  variables   in  YL(l). 
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ERROR(l)    for  I   =  l(l )NY  stores   the  sum  of  the   corrections 
computed  during  each  integration  step.      It   is   initialized  to   zero   after 
the  predictor  evaluation   and  is   incremented  by  Fl(l) — the  improved 
corrector — after  each   corrector  step.      If  the  integration   step  is 
successful,    all  elements   of  Y(J,l)    are   updated  by  Y(j,l)   =  Y(J,l)    + 
A(J)*ERR0R(I).      The    difference    for  ERROR(l)   between  two  integration  steps 
is   used  to    calculate   the  next  higher  derivative   element  when  the  best 
step  size    for  order  k+1  is  being  determined.      When  it   is  known   that  on 
the  next   step   a  new  value   of  H  will  be   calculated,  ERROR(l)    is   saved  in 
ERSV(I)    for  I   =   1(1)NY. 

YMAX(I)    for  I   =   l(l)NY   stores  the  maximum  of  1.0  or  the 
largest  value  of   |Y(l,l)|    computed.      YMAX(l)    is   used  in  the  error 
evaluation   for   an   absolute  error  test  when  Y(l,l)   has  never  exceeded 
i.O   and  a  relative  error  test  when   it   has. 

A(j)  ,  J  =   l(l)NQ+l,   stores   the  updating  factors   for  the  higher 
order  derivative   terms  in  Y.      It   corresponds    to  the  vector  2_  in  Section  2, 

G(J)    is   storage   for  global  variables,   parameters   that  may  be 
easily   changed  from  one   integration  to  the  next.      They   allow  a  user  to 
experiment  with   different  parameters   in   attempts   to  optimize    given 
systems   under  simulation. 

T(J)    contains   variables   that   are   dependent  only  on   G(*)    and 
the  independent  variable  T(l).      Thus,   they  need  to  be   evaluated  only 
for  each  value   of  the   independent  variable. 

PW(I,J)    is   the   matrix  J. 

The  remaining  calling  arguments  to  DFASUB  follow: 


12 

EPS  is  the  error  per  step  criterion  that  nust  be  met  by  the 
L2  norm  of  the  error. 

EQN  and  VAR  are  two  arrays  that  are  required  by  routine 
MATSET  (see  following)  when  MATSET  is  compiled  by  the  general-purpose 
system  [l].   When  this  is  not  the  case,  they  can  be  dummy  arrays. 

HMAX  is  the  largest  step  size  that  DFASUB  will  take.   It 
should  be  chosen  less  than  the  period  of  any  known  periodic  solution. 
Otherwise,  no  special  care  need  be  taken  in  choosing  it. 

HMEN  is  the  minimum  step  size  taken.   H  is  the  initial  step 
size  used  unless  it  causes  convergence  problems  ,  in  which  case  a  new 
step  size  is  chosen. 

JSTART  is  a  start  indicator.   When  it  is  zero  on  a  call  to 
DFASUB,  the  routine  initializes  itself.   When  it  is  1,  it  assumes  that 
initialization  has  taken  place  and  continues  integration.   On  exit, 
it  contains  the  current  order. 

KFLAG  is  the  output  completion  code.   Its  meanings  are 

+1  successful  integration 

-1  error  test  failure  for  H  >  HMEN 

-2  too  many  floating  point  errors 

-3  the  corrector  failed  to  converge  for  H  >  HMIN 

-k     the  corrector  failed  for  even  the  first  order  method 

It  is  possible  to  write  into  the  main  program  a  routine  to 
try  possible  remedies  when  DFASUB  quits  with  an  unsuccessful  completion 
code  making  use  of  KFLAG  and  the  computed  GO  TO  in  FORTRAN.   After  these 
remedial  measures  are  taken,  JSTART  should  be  set  to  1  and  DFASUB  called  again. 
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M  is   the   dimension   of  DY ,  normally  the  same   as  N. 

MAXDER  is  the  maximum  order  method  used,   never  more  than  6. 

ML   is   the  number  of  variables   in  Y(l,J),   J  =   l(l)ML,  which   are 
included  in  the  error  test.      If  the  error  in   any  group  of  variables   is 
not    critical,  then  they  should  be  placed  at  the  end  of  the  Y   array   and 
ML  set  to  exclude   them  in  the  error  test.      N  is   the  total  number  of 
equations.      NL  is  the  number  of  linear  equations. 

TEND   is  the    final  value   of  T(l).      Integration  stops  when  this 
value   is   reached. 

Subroutines 


The   dependent  variables   in  DFASUB   fall   into  three    classes  : 
those   dependent   only  on   global  terms  or  parameters   in  the   G(*)    array, 
variables    described  by  linear  equations  which  have  Jacobian  elements 
dependent  at  most  on  the   arrays    G(*)    and  T(*)    so  these   are   the   elements 
of  the  matrix  P  of  Section  2,    and  variables  whose  entries   to  the 
Jacobian  matrix  change  whenever  some   other  variable    changes   its   value. 
To  handle  these    differences   efficiently,    five   subroutines   called  by 
DFASUB   are    concerned  with   computing  the  Jacobian  matrix  and  performing 
the   corrector  iteration    (2.3). 

Since  the  elements   in  the  Jacobian  have   different   dependencies, 
they   can  be   recomputed  at  different  times  by 
MATSET    (A(2)  ,DY ,EPS ,EQN ,G,HINV,0,M,MF,N ,NY ,IT,PW,F1 ,T ,VAR,Y,YL) . 

When   IT  =   1,   only  those   values   in  the  Jacobian  that   are   fixed 
for  one  set  of  parameters   G(*)    are    calculated.      This   need  be   done  only 
once  per   integration.      When   IT  =   2 ,   only  those  values   that   are   dependent 
on  T(*)    and  G(*)    are   computed.      This   is    done  whenever  a  new  step   is 


Ik 

started,    causing  a  new  value   T(l)  =  T(l)  +  H  to  be    computed,   and 

3  new  O-Lu. 

whenever  the  step   size   is   reduced  and  the  integration  step   restarted. 
When  IT  =   3 ,   all  remaining  elements    are    computed.      This  must  be   done 
at  each   corrector  iteration. 

The   corrector  iteration    (2.3)    does  not    actually  multiply  by  J 
but  rather  performs   a  symbolic  inversion  of  the   sparse  matrix  PW.      Again, 
the  symbolic  inversion   of  PW  is  broken  up   according  to  the   dependence  of 
each  element   on  T ,    G,    and  Y.      MATINl(PW)    inverts   just  those   elements 
that   are   dependent   only  on   G.      MATIN2(PW)    inverts   elements    dependent  on 
T  and  G.      MATIN3  inverts   all  remaining  elements. 

MATMUL(PW,DY,Fl)   performs   the  equivalent  of  solving 
(PW)(Fl)   =  DY   for  Fl;   DY  holds   the  values   of  F(y,   v)   =   F(y ,  y' ,   t)    +  P(t)v; 
and  Fl  is   a  vector  used  in  the   correction  step. 

DIFFUN(T,G,DY,Y,YL,HINV)    is  the   differentiating  subroutine  which 
places  y!    -   f .  (y_)   =  HINV*Y(2,l)    -  f(y(l,l))    into  DY(l)    for  I  =  l(l)NY, 
and  the   correction  to  the   linear  equations   into  DY(l),   I  =  NY+l(l)N. 

The  remaining  subroutines    are  not   connected  with  the 
differential  algebraic  equation   system,  but   are   used  for   convenience. 

KNTSPl(l)    is    an  assembly  language   function  that  keeps   track  of 
floating  point  overflow  and  underflow,  and  returns    a  value   of  three   or 
more   if  enough  such  errors    are   encountered  to  invalidate   the  results   of 
th  e   co  mput  at  i  on . 

COPYZ(S,X,L)    copies  the  L  values   of  single  precision  array  X 
into  single  precision   array  S.      The  actual  variables   used  are   double 
precision   arrays   so  L  is   twice   the  size   of  the   actual  arrays.      If 
X  =  Y,   S  =  SAVE,  we  have   L  =   lU*NY.      For  X  =  YL ,   S  =  YLSV,  we  have  L  =  2*NY. 

S2(T,G)    evaluates   those  variables   stored  in  T(*)   that   are 
dependent  on  T(l),   the   independent  variable;   and  G(*),  the    global  parameters 
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h.      USER  PROGRAM  FOR  DFASUB 

To  write   a  program  to  use   DFASUB   all  the   subroutines    called 
by  it   must  be  present  in  some   form.      The   subroutine  DIFFUN ,  which  is 
supplied  by  a  system  compiler  in  the   general -purpose   system,   must  be 
provided  to  solve  the   system 

f(ys  y1,  t)  +  P(t)y_  =  0. 
The   calling  arguments    are   the   T(*)    and  G(*)    arrays,    the  Y   array    containing 
y.    in  Y(1,I)    and  H*y!    in  Y(2,l),   HINV  =   l/H,    and  v.    in  the   array  YL.      The 
system  can  be  put  into  the  proper  form  by  expressing  the   equations    as 

0  =  yl  -(f1(Y»xL,T,G)   +  P^) 

0  ■  yNY   "( VY,YL,T,G)    t   ?mZ) 

°  7    W(Y'YL'T'G)    +  PNY+1^ 
0   =    fN(Y,YL,T,G)    +   PNv 

for  P     the   I-th  row  of  P(t). 

The   zeros   are   then   replaced  by 

DY(I)  ,   I   =   1(1)N. 
An  example  system  is   to  be  found  in    [5].      It  is 


0  =  y!    -  s   +    (r-y. )      +     E     b. .   y.,  i   =    l(l)U 


i=l 


ij      J 


for 


r  =  (y1  +  y2  +  y3  +  yu)/2 


s  =     E      (r-y    )"/2 
i=l 


[b      ]-B 


i+UT.5+e  -1+52. 5+e  -^7.5+e  -52.5-e 

■i+52.5+e  UUT.5+G  52.5+e  U7.5-e 

-HT.5+e  52.5+e  UU7.5+e  ^52. 5+e 

-52.5-e  UT.5-e  ^52 . 5-e  UUT-5+e 


for   e   =    .00025. 
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0  =  *5  +  yl  y6  +  yl  y6 

0  =  2H  =  y6  "  yi  +  vi  "  x  "  e_t 


=  vl  +  F6 
0  =  v1  -  v2  +  yx  y6 

=  vx  -  v2  +   FT 

0  =   v     +  v     +   5v     v 
1  2        ^1  y2 

=  v1  +  v2  +  F8 

for  y<°>    =  y<°>   -   1,   v[°)   =  -2,   v^   =   -3,  y[0)   =  -1,   i  =   l(l)U. 

Appendix  I    contains   the   FORTRAN   routine   DIFFUW(T,G,DY,Y,YL,HINV) 
that  is  used  to  evaluate   this    system  as  veil  as   a  sample   routine  S2(T,G) 
to  solve   T(2)    =  EXP(-T(l)). 

The   computation   of  the  Jacobian   can  be  split  up  into  the  three 
subsets,   each  of  which  is   evaluated  by  MATSET   depending  on  the  parameter 
IT.      Three   different  inversions   are   also   required.      If  the  user  would 
rather  not  write  these   four  routines,    a  simpler   routine  MATSET   can  be 
written  so  that  nothing  is    done  when  IT  is   equal  to   1  or  2  ,    and  the  entire 
Jacobian  matrix  is   generated  when  MATSET   (...,3,...)    is    called.      The 

evaluation   of  the  Jacobian   can  be    done  by   a  special  routine   that   computes 

9f.  9f. 

- —  and  - — ; —  according  to   an   analytic   formula,    or  else  numerical 
yj  y  J 

differencing  as   in    [6]   requiring  N   calls   to  DIFFUN  may  be   used.      MATIN 3 
can  be   replaced  by   a  call  to  a  standard  Gaussian  elimination  subroutine, 
while  MATIN  1  and  MATIN2   are  made   dummy  subroutines.      MATMUL   can  be   replaced 
by  a  back  substitution   subroutine  written  to  work  with  the  routine   in  MATIN3. 


IT 

Two  implementations   of  this  procedure   are  possible.      The 
easiest   one  is  to  write   a  subroutine  which  has   a  place   in  its    calling 
sequence   for  each  of  the  parameters   in  the  calling  sequence   in  DFASUB , 
and  then  manipulate  these   parameters,  possibly  preparing  them  for   a  call 
to   a  second  subroutine.      This  has  been   done   in  the  remainder  of 
Appendix  I.      For  example,   MATSET  has    dummy  variables    for  EQN   and  VAR , 
two   arrays   MATSET   as    compiled  by  the   general-purpose   system  use  but   are 
not  needed  in  this   version,  which  performs   numerical   differencing  to 
compute  PW  when  IT  is    3  and  does   nothing  when  IT   is   not    3.      MATIN 3 (PW) 
simply   calls   Moler's   routine    [7]   DECOMP.      Since   DECOMP   requires   the 
number  of  variables  N,   this  has   been  provided  in  the  unnamed  COMMON  block 
in  MATMUL   and  DFASUB.      MATMUL   calls   Moler's   routine   SOLVE  which   also 
needs   N.      Both  of  these   routines    communicate  through   an   array  IP,   which 
is  provided  in  a  named  COMMON/IPP/    and  dimensioned  at   least   as   large 
as   the  system. 

The   second  implementation  requires  that  the  user  remove   all 

superfluous   subroutine    calls   and  change  the    calling  sequence   of  the 
routines  he   does   use,  thus   saving  the  overhead  of  dummy  subroutine   calls 
and  extraneous   preparation   of  data.      This  has   the   disadvantage  that   if 
DFASUB   is   updated,    all  this   work  must  be   repeated  on   the  new  version; 
whereas   if  the  subroutine    calls    are   left  alone,    those   that  work  with 
one  version   of  DFASUB  will   probably  work  with   any  new  version. 

Usable  versions   of  COPYZ   and  KNTSPI    are   also   included  in 
Appendix  I.      If  the   computer  in   use  has   some  way  to   count   suppressed 
underflow   and  overflow  under  program  control,   KNTSPI    can  be   programmed 
to   return   a  value   greater  than    3  when  too  many  errors   occur. 


18 

Finally,    DFASUB   must  provide   the   initial  values    of  Y(l,l), 
Y(2,l),   and  YL(J).      The   general -purpose   system  calls   a  program  DIFMF3    [8] 
which,   if  given  either  Y(l,l)    or  Y(2,l)    for  each  I   and  YL(j)  ,  will   find 
all  necessary  initial  values.      If  all  initial  values   Y(l,l)    are  known , 
then  by  setting  Y(2,l)    =  0   for  I   =   1,NY,    calling  DIFFIM  once   and  setting 

Y(2,I)   =    -DY(I), 
good  approximations   can  be   given  to  all  necessary  initial  values. 
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5.      PROGRAM  LOGIC 

An   anthropomorphic  view  of  the   actions   of  DIFSUB  while  solving 
a  system  of  equations    follows.      A  main  program  specifically  written   for 
the  PDP -8/ GRAPHICS   system   [l]  provides   the   calling  sequence   that  enters   the 
subroutine  DIFSUB  and  the   subroutines  which   DIFSUB   calls   are   compiled  by 
the  system  loaded  with  DIFSUB.     A  detailed  description  of  what  each   of 
these   do  is  to  he   found  in  Section   3. 

The  program  is   entered  with  JSTART  =  0.      The  bookkeeping  variables 
that   record  the  number  of  steps,  NS ,    and  the  number  of  matrix  evaluations,  NW, 
are  initialized  to  0.      Both  y  and  y'    are   expected  to  be   in   array  Y(l,l)    and 
Y(2,l).      This    can  be    accomplished  by   a  call  to  the  steady  state  problem 
package   DIFMF3    [8].      Y(2,l)    is  scaled  by  H,  YMAX(l)    is   set   to  1.0,   the   order 
NQ  is   set  to   1,    and  the  variables  l_.    are  placed  in  A(j)  ,   J  =   l(l)NQ+l. 

J 

Various   error  test   criteria  are   computed  in  E,   EUP  ,   and  EDWN ;   ENQ1 ,  ENQ2  , 
and  ENQ3  held  l/(2*NQ),   l/(2*(NQ+l),    l/(2*(NQ+2)    the  power  to  which  the 
ratio   of  computed  error  to  desired  error  must  be   taken   to   find  the  ratio 
of  new  to   old  step  size. 

MATSET  is   called  with   IT   =  1  to  evaluate   any   constant   elements   of 
PW;   MATIN1  inverts   these  elements. 

The   value   of  Y(J,l),   J   =    l(l)NQ+l,   I    =   l(l)NY    are  saved  in 
SAVE(J,l)    in   case   they  are  needed  to   restart  the  program;   YL(l),   I   =   l(l)NL 
is   saved  in  YLSV(l).      H,   T,    and  NQ  are  also   saved.      DIFSUB   then    calls   S2 
which  evaluates   any  variables   there   are    functions   only   of  T   and     G.      MATSET 
is    called  with   an   argument  of  2   to  evaluate   and  update   any  part   of  the 
Jacobian   matrix  PW  which  would  change   as    a  result  of  the   call  to  S2.      DIFSUB 
then   calls   MATIN2  which   does    an   inversion   of  parts   of  PW  that    change    as    a 
result  of  calling  S2   and  MATSET    ( .  .  .  ,2  , . .  . ) . 
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DIFSUB   multiplies    Y(J,l)   by  the   Pascal  triangle   matrix  as 
described  in    [6].      This    replaces  y.        with  y.     ,    s,   the  predicted  value 
of  the  nonlinear  variables. 

The   corrector  step  is   a  long  loop  that   starts  by  initializing 
ERROR   to  0.      DIFSUB    calls   DIFFUN (T, G,DY ,Y ,YL,HINV)   which  places    in   DY 
the  updated  y!    ,    v   value.      MATSET    (...,3,...)    is  now   called  to  update 
the  Jacobian   as    a  result   of  having  the  predicted  value  of  y.    and  the 
corrected  value   of  y!    available.      MATIN3(PW)    is    called  to   make   any 
changes   in  the  LU  decomposition    as    a  result   of  the   call  to   MATSET 
(...,3,...).      IWEVAL  is   set  to  -1  so  that   all  of  this   re-evaluation   of 
PW  need  not  be   done   again  unless   it   seems   necessary.      MATMUL(PW,DY ,Fl) 
is    called  which   does   the  back  substitution   that   computes   the    amount  by 
which  Y(j,l)    and  YL(l)    are  to  be    changed  to  hold  the   corrected  value. 
YL(I),Y(1,I)    and  Y(2-,I)    are    corrected  by  A(J)*Fl(l )  ;    ||Ay,    J|      is 
then   computed.      If   |  |Ay|  |      is    less   than  BND,   the   computed  error  bound 
for  the   order,   then  the  program  continues.      If  not,   MATMUL   is    again 
called,  YL(l),   Y(l,l),   Y(2,l)    are    corrected,  Ay ,    \    is   added  to  Ay/-,  \ 
and  placed  in  ERR0R(*). 

Actually,   not   only  the   present   corrector  term  but   also   an 
estimate  of  the  next   corrector  term  is    compared  to  BND.      Since   the  L2 
norm  of  the   i-th   corrector  term   |  |  ERROR    |  |      is    approximately 
R*|  |  ERROR        I  |      for  R   constant  throughout   each   small   region  of 
integration,   the  predicted  next  value   of  the   corrector  can   also  be 
compared  to  BND   and  the  method  is    considered  to  have    converged  if 

MIN(| | ERROR1 | | ,    R*| |ERR0R1+1| | *2)    <BND. 
R  is   initially   set   to    |  |  ERROR    |  \/\  |  ERROR    |  |    and  updated  at   each  step. 
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Assume  that  the  corrector  term  did  converge  in  two  corrections. 
Then  the  estimated  error 

|  |error(i)/ymax(i)|  I 

is    computed;    and  if  this   is   less   than  E,   the  maximum  error  hound  allowed 
to  have  the  real  error  per  step   less   than  EPS,   the  program  continues. 
Otherwise,    a  better  value   of  H   for  this   or   one   lower  order  is    computed 
as    described  below,  KFLAG  is   reduced  by  2,   the  original  values   for  this 
step  are   returned,    and  the   step  is    repeated  for  the  new  step  size.      Should 
KFLAG  reach  -5,   indicating  three  steps   taken  with  three   different   step 
sizes  yet  none   gave   an   acceptably  small  error,    the   order  is   reduced  to   1 
and   another  attempt   is   made.      If  this   is    unsuccessful,    KFLAG  is   set   to 
-k  and  DIFSUB  makes    an  error  return   to  its    calling  program.      If  at   some   time 
H  becomes    less  than  HMIN,   KFLAG  is   set  to  -1  and  DIFSUB  makes    an   error  return 

Otherwise,  the  rest  of  the  Y(J,l)  array,  if  NQ  >  1,  is  updated 
using  the  contents  of  ERROR(l)  by  Y(J,l)  =  Y(J,l)  +  ERROR( I ) *A( J) ,  then 
DIFSUB  returns  to  the  beginning  of  the  program.  JSTART  has  been  set  to 
1  so  no  initialization  occurs.  The  integration  step  procedes  as  in  the 
first  step,  except  that  after  NQ  successful  steps  ERROR(l)  is  stored  in 
ERSV(l)    and  the  next   step   (if  successful)    uses   this    and  other  information 

to   check   for   a  better  step   size   and  order. 

(k) 
Since  we   also  have   an  estimate   for  y      and  y        /k!    stored  in   the 

n  n 

array  Y(7,l),   and  since  ERROR(l)    stores   the   current  step's    accumulated 

corrections    for  Y(l,l),   the  backward  difference   of  the   last    component   of 

Y(K,N)    is  ERROR(l)*A(K)   which   is    an   estimate    for 

,k+l      (k+l) 
h         y 
_n 

k! 


22 


Also,   since  ERROR(l)    corresponding  to  the   last  y  are  saved  in  ERSV(l) , 

we  know  that 

ERSV(I)    -  ERROR(I)    =  h(k+2)   y(k+2 } (tQ) /k! . 
These  values   can  he  used  to    compute   the  "best   step  size   and  order  in 
the  following  ways . 

After  the  predictor-corrector  method  converges,  we   can   check 
to   see   if  the  truncation  error  was  within  hounds   by  seeing  if 

2          Ml     Ck+1  b^  y<^>    g 
EPS     1  £   '  YMAX(I) > 

Ml      C.       *ERROR(l)*K!*A(K) 

=      y      f    k+1 \^ 

*      v  YMAX(I)  y 


by    comparing 


to 


E  =    (EPS/C.  ._    A(K)*K!)2 


Ml 
_  v^    ,ERR0R(lK2 

D  "  zz    4max(i)   j   ' 

When  we  want  to  see   about    changing  order  and  step   size    (after 
K+N*9,  N  =  0,1,...,   successful  steps   at   order  K   and  step  size  H)  ,  we   can 
compute  the  best   step  size   to  use   at  the   current   order,   order  1  higher 
and  order  1  lower.      Since 

D  ^    fHNEWs2    (p+1) 
E     '    [     H     ' 

by   computing  PR2   =   l^D/E)1'2^"*"1' , 

MEW  =   H/PR2 

will  give   the  best  step  size   for  order  K   (1.2   is     a  heuristic   constant 

k+2 
to   compensate   for  ignoring  the  0(h        )   terms). 
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letting 


and 


If  one   order  lower  method  is   used,  then 

C      hk   v'k^  C    *v(k   l)*r>' 

EPS     ^Z(      YMAX(l)    }      =E(        YMAX(l)       } 


EDWN   =    (EPS/C   *k! )2 


then 


D  "  MYMAX(l)J 


PR1  =   1.3(D/EDWN)1/2p 


gives  the  factor   for  the  "best  step  size   for  order  k-1  methods 
Finally,   if  order   one  higher  is  used,  we  need 

n         >k      (k) 
2  Ck+2  h     yI      s2 

^    ^      mx(i)     } 

C.     Q    (ERSV(I)-ERR0R(I))A(K)*K!) 
v  YMAX(I) 

Letting  EUP  =    (EPS/C        *A(K)*p!)2  and 


k+2 


then 


and 


_     .    wD/ERR0R(lh2 
D  "    M      YMAX(I)    j 


D  _    ,HNEW,2    (p+2) 
E   "     ^    H     ' 


pes  =  i.m^)1/2  (p+2) 

gives   the  best   step  size.      The   1.3  and  l.k  terms   bias   the  method  toward 
not   changing  the  order  or  picking  one  order   lower,    respectively,   since 
these  would  minimize   overhead  calculations. 
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Once   a  new  step  size   has   been   computed,   the   array  Y(J,l)    must 
be   changed  to   reflect  the  new  step   size   since  Y(J,l)   =  Y  h  /(J-l)!. 

This    is   accomplished  "by  multiplying  the  elements    of  SAVE ( J, I )    (if  the 
single   step  error  were  too  large   and  the  step  is  being  repeated)    or  of 
Y(J,l)    (if  a  new  step  size   is  being  prepared  for  the  next  step)   by 
(H/HOLD)(J_l). 

This   is    done  in   a  DO  loop   after  statement   750   or   800   depending 
on  whether  SAVE  or  Y   is  supplying  the  values   to  be  reached.      If  no 
significant  improvement  in  step  size   can  be  made,   IDOUB   is   set  to   9   (or 
other  large  integer)    and  the   step  size   is   re-evaluated  again  if  the 
steps    continue  to  be  successful. 

If,    contrary  to  our  assumption,   the    corrector  steps    do  not 
converge  after  two  tries,   then  H  is   reduced  to  E/k  unless   IWEVAL  =   0 
indicating  that  PW  should  be   re-evaluated  first.      Y(J,l)    is   scaled  to 
reflect  this  new  step  size    and  a  new  integration   step   is   attempted.      If 
H   cannot  be  reduced,  KFLAG  is  set   to   -3   and  DIFSUB  makes    an  error  return. 

Each  time    a  new   step  is   entered,   DIFSUB   checks    for  T  >   TEND 
and  KFLAG  <   0.      If  either  of  these   is   true,   the  program  exits   to  the 
calling  routine  with  JSTART  set   to  the  present  value   of  the  order  NQ , 
HNEW  set  to  the  best  H   for  the  next   step,   H  set  to  the  present  H,    and 
Y(j,l)    set  to  its   last   successful  value. 

The  entry  point  REDSUB   is  used  to   re-enter  the  program  after 
problems  with  the  matrix  processing  routines. 

This   section,  while   describing  the   general  program  logic,    leaves 
out  many  details   of  programming  that  would  be   unwieldly  to   describe  here. 
The  theory  and  some   ideas    about  the  programming  techniques    are   in  the   other 
sections,    and  after  reading  these   and  using  the   program  copy  in  Appendix  II, 
this   section  should  be   sufficient   to  understand  the  program  logic. 
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SAMPLE  PROGRAM 


27 


c**** 

c* 

c**** 


c**** 

c* 

c**** 


10 


20 


c**** 

c* 

c**** 


30 
35 


999 
50 


IMPL 
DIME 

1  G(l 
DATA 

«•  -45 

♦  52. 

♦  452 
***** 

SUPP 

***** 

DATA 

♦  5.D 
CALL 
CALL 

***** 
SET 

***** 
nr  i 

Yd, 
Yd 

Yd, 

DO  2 
Y(2, 

YL(  1 
YLI2 

M  =  N 
***** 
FIND 
***** 
CALL 
DO  3 
Y(2, 
JSTA 
CALL 

1  JST 

2  YL 
WRIT 
FORM 
STOP 
FND 


AL*8( A-H.Q-Z) 

(7,6),YL(2),SAVE(7,8),YLSV(2),PW(R4), 

, 0Y(8) , ERSV(6) .ERROR (6) ,F1( 8) ,EQN(1 t ,VAP  (1) ,YMAX( R) 

5002  5,-452.49975,-47.499  75,-52.5002  5, 

,44  7.50025,52.50025,4  7.49975,-47.49975, 

47.5002  5,4  52.499  75,-52.5002  5,47.49975, 

447.50025/ 

********************************* 

ERFLOW  AND  UNDERFLOW 
********************************* 

,EPS,HMAX,HMlN,H,T,TENP/8,2,8,l.D-4, 
0,  l.D-4,0. rO, L. DO, 1.03/ 
(208,256,-1,1 ) 
(207,256,-1,1 ) 

********************************* 

VALUES    AND    ZERO    FIRST    DERIVATIVES. 
********************************* 


ICIT  RE 
NSIPN  Y 
6),T(2) 
G/447. 
2.49975 
50025,4 
.49975, 
******* 

KESS    OV 
******* 

N,NL,M 
2,1.D-1 

ERPSET 

EPRSET 

***  **** 

INITIAL 
******* 

0     1=1,4 
I)=-l. 
,5)=1. 
6)  =  1. 
0    1=1,6 
I »=0.D0 
)=-2. 
)=-3. 


**************************************** 

FIRST    DERIVATIVES    OF    N1NLINEAP    VARIABLES 
**************************************** 

DIFFUN(T,G,DY,Y, YL , L.I 
0    1=1,6 
I )=-PY( I ) 
RT  =  0 

DFASUB(DY,FPS,EQN,ERPOR,FRSV,Fl ,G, H, HMAX, HMIN , 
ART,KFLAG,M,6,6,N,NL,PW,SAVE,T,TEND,VAR,Y,YL, 
SV, YMAX) 
E(6,999)KFLAG 
ATI  '      KFLAG  =•  ,16) 


SUBROUTINE  S2IT.G) 
REAL*8  T«2) 
T(2)«DEXP(-T(1I I 
RETURN 
ENO 


SUBROUTINE  DIFFUNJ T, G ,DY, Y, YL , HI NV ) 

IMPLICIT  REAL*8  (A-H.Q-Z) 

DIMENSION  G(4,4),DYI8),Y(7,6),YL(2>  ,T(2> 

R=»(Y(l.l)*Y(l,2)*Y(l,2>*Y(l,4)l/2. 


S  =  0. 

DO  20  !«lt* 
20  S  =  S*(R-Y(l, I ) »  **2/2- 

DO  3  0  1=1,4 

DY( I l=HlNV*Yt2,II-S*(R-Y( 1,  I)  )*»2 

00  2  5  J =1,4 
25  DY< I  )=DY(I  )*G( I , J  > *Y (  1, J) 
30  CONTINUE 

0Y«5)=HINV*(Y( 2,5)*Y(  1, ll*Y( 2,6)*Y(2,l)*Yl  1,61  I 

DY(6)*2.*Y(  l,6)*Yd,6)»*3-Y(  I,  ll*Ylt  1I-1.-TC2I 

0Y<7)=YL( 1)-YL(2)*Y(  l,l)*Y( 1,6) 

0Y(8)  =  YLd  )*YL  (2)»5.*Y(1,1)*Y(  1,2) 

RETURN 

END 
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SUBROUTINE  CDPYZ(S,Y,L) 
DIMENSION  Y(  1)  ,S (  1  ) 
DO  10  J  =  1,L 
S(J)=Y(J) 

RETURN 
END 


FUNCTION  KNTSPim 

KNTSPI=I 

RETURN 

END 


SUBROUTINE 

RETURN 

END 


MATINKP) 


SUBROUTINE  MATIN2IPI 
RETURN  'IK| 

END 


C* 


C* 


c* 


c* 
c* 
c* 


c* 


c  + 


SUBROUTINE  M AT  SET ( A2 , CY , EPS , D4, G,HI NV, D5 ,06 , D7, Nt NY, I CHK , 
1   PW,F1,T,D9,Y,YL) 

IMPLICIT  REAL*8( A-H.Q-Z) 

DIMENSION  OYU), G(l),  PW  11), TU). Y(7,1),YL(1),F1(1) 

SEE  IF  THIS  IS  A  CALL  TO  BE  FULLY  HANDLED. 

IFI  ICHK.NE.3)  GO  TO  100 

NL=N-NY 

NN=N*N 

DO  20  1=1, NN 
20  PW( I )  =  0. 

DO  40  J=1,NY 

SAVE  COLUMN  ELEMENTS 

F=Y( 1, Jl 

E=YI2,J) 

R=FPS*DMAX1 (EPS, DABS  (  F)  ,DABSIG) ) 

FIND    A    DIFFERENCE    ELEMENT    =    MAX( EPS**2 , EPS*Y( 1 , J ) , EPS*Y ( 2 , J ) ) 

Yd,  J)  =  Y(1,  J|«-R 

Y(2, J)=Y<2, J)-A2*R 

CALL    DIFFUN(T,G,F1,Y,YL,HINV) 
DF  A(2)     DF 

ASSIGN  VALUES  TO  --  

DY      H   DY» 

DO  30  I=1,N 
30  PW< I*(J-l)*N)= IFK I)-DY( I»)/R 

RESTORE  COLUMN  ELEMENTS 

Y(2,J)=E 
40  Y< 1, J |=F 

IF  ANY  LINEAR  ELEMENTS,  FIND  DF/DYL. 

IF(NL.E0.0)GO  TO  100 

DO  80  J=1,NL 

F=YL( J) 

R*EPS*0MAX1IEPS,DABSI F) » 

YL( J)=YL(J)*R 

CALL    DIFFUN(T,G,F1 ,Y,YL,HINV) 

DO  70  1  =  1, N 
70  PW(  I*JJ*NY-1)*N)»IFH  II-DYI  I)  l/R 
80  YL<J)*F 
100  RETURN 

END 


SUBROUTINE    MATIN31PW) 

IMPLICIT    REAL*8     (A-H.Q-Z  » 

DIMENSION  PW(1 ) 

COMMON  N 

COMMON     /IPP/    IP(40» 

CALL    DECOMP(N,N,PW,IP) 

RETURN 

END 


in 


SUBROUTINE    M AT MUL( PW , DY , Fll 

IMPLICIT    REAL*8     (A-H.O-Z) 

COMMON    N 

COMMON  /IPP/  IPUO) 

DIMENSION  PWI1),DY(1)  ,FIC 1) 

DO  10  I-1,N 

Fl(  I  )*DY(I) 

CALL  S0LVE(N,N,PW,F1, IP) 

RETURN 

END 
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SUBROUTINE  CECOMPI N , NCI M,  A  ,  IP ) 

IMPLICIT  REAL*8<R-H,Q-Z) 
C 

C     MATRIX  TR IANGULARIZATICN  BY  GAUSSIAN  ELIMINATION. 
C 

C    INPUT... 

C  N    =    ORDFP    OF    MATRIX 

C  NDIM    =    DECLARER    DIMENSION    OF    ARRAY    A. 

C  A    =    MATRIX    TO    RF    TR I  ANGULAR  I  ZED.     (FOR    STIFF    METHODS,     A     IS    SINGLE 

C  PRECISION;       ALL    OTHER    VARIABLE    AP F    OOUBLE    PRECISION.) 

C         OUTPUT... 

C  A(I,J),     I.LE.J       =    UPPER    TRIANGULAR     FACTOR,     U. 

C  All, J),  I.GT.J  =    MULTIPLIERS       =    LOWER    TRIANGULAR    FACTOR,     I-L. 

C  IP(K),       K.LT.N    =    INDEX    OF    K-TH    PIVOT    ROW. 

C  IPIN)       =    «-l)**INUMBEP    OF    INTERCHANGES)    OR    0. 

C         USF     'SOLVE'     TO    OBTAIN    SCLUTION    OF    LINEAR     SYSTEM. 
C         DPTERM(A)    =     IP(N)*A(1,1)*A(2,2)*...*AIN,N). 
C  IF     IPIN)    =    0,     A    IS    SINGULAR,    SOLVE    WILL    DIVIDE    BY    ZERO. 

C 

DIMENSION    AINDIM.NDIM), IP(NDIM) 

DIMENSION    BI4,*) 

COMMON    NFNS,NW,B 

NW    =    NW    +    1 

IPIN)    =    1 

DO    6       K=1,N 

IFIK.FQ.N)    GO    TO    5 

KP1    *    K+l 

M    =    K 

DO    1        I=KP1,N 

IF(ABS(A(I,K)).GT.ABS(A(M,K)I)  MM 

1  CONTINUE 
IP(K)     =     M 

IFIM.NE.K)        IPIN)    =    -IPIN) 

T    =    A(M.K) 

A(M,K)     =    AIK.KI 

AIK.K)    *    T 

IFIT.EQ.O)    GO    TO    5 

DO    2       I*KP1,N 

2  A( I,K)    =    -All ,K)     /    T 

?°    \,i!n1,N  SUBROUTINE     SOL  VE  IN.ND  IM  ,  A,  B,  IP  ) 

AIM,;!    ^AIK.J,  c  ,MPUCIT    REAL*8     'B-H.O-II 

Jf^EoIoI)    GO    TO    4  I  S°LUTI0N    °F    L,N6AB    SYSTEM'    A*X    '    B« 

DO    3         I-KPUN  c  ,NPUT... 

I  rn^JM,^'    "    AII'JI    *    ACI'K,*T  C  N    -    ORDER   OF    MATRIX. 

*  CONTINUE  c  N0,M    ,    DECLARED    DIMENSION    OF    ARRAY    A. 

5  IF(A(K,K).EQ.O.)     IPIN)    -0  C  A    «    TR I ANGUL ARI ZED    MATRIX    OBTAINED    FROM    'DECOMP' 

6  CONTINUE  C  B    ■    RIGHT    HAND    SIDE    VECTOR. 

RETURN  C  IP    -    PIVOT    VECTOR    OBTAINED    FROM    'DECOMP'. 

END  C         OUTPUT... 

C  B    ■    SOLUTION    VECTOR,     X. 

DIMENSION    AINDIM.NDIM),     BINDIMI,     IP(NDIM) 
IFIN.EO.l)    GO    TO    9 
NM1    -    N-l 
DO    7      K"1,NM1 
KP1    ■    K    ♦    I 
M    ■     IPIK) 
T-BIM) 
BIM)     -    P(K) 
BIK)    -    T 
DO    7       I-KP1.N 

7  BID    -    B(I)    ♦    AII,K)*T 
DO    8         KB-l.NMI 

KM1    -    N    -   KB 

K    ■    KMl    ♦    I 

BIK)     ■    BIKJ/AIK.K) 

T    *    -BIK) 

DO    8       I-1.KM1 

8  Bl  I  )    «    BID    ♦    A(I,K)*T 

9  BID    -    BU)/AI1, 1) 
RETURN 

END 
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SUBROUTINE    OFASUB    (  DY  ,  EPS  ,  EON,  ERP  OR  ,  EP  S  V, 
+  Fl ,G,H,HMAX,HMIN, 

♦         JSTAPT,KFLAG,M,MAXDER,M1 ,N, 
«■         NL, PW, SAVE, T, TEND, VAR,Y, 
«■  YL.YLSV, YMAX) 

IMPLICIT    PEAL**    (A-H.O-Z) 

C* 

(-*     THE  PARAMETERS  to  THE  SUBROUTINE  PIFSUB  HAVE 

C*     THE  FOLLOWING  MEANINGS: 

C* 

C*  N         THE  NUMBER  TF  VARIABLFS. 

C*  NL       THE  NUMBER  OF  LINEAR  VARIABLES 

C*  (NY  =  N-NL   IS  THE  NUMBER  OF  VARIABLES  WITH  DERIVATIVES) 

C*  Ml       THE  NUMBER  OF  EQUATIONS  TO  TAKE  PART  IN  THE  ERROR  TEST, 

C*  TFND     END  CRITERION 

C*  T         THE  INDEPENDENT  VARIABLF. 

C*  G        AN  ARRAY  OF  GLOBAL  VARIARLES 

C*  Y        A  7  BY  NY  ARRAY  CONTAINING  THE  DEPENDENT  VARIABLFS 

C*  AND  THEIR  SCALED  DERIVATIVES.  Y(J*l,I)  CONTAINS 

C*  THE  J-TH  DERIVATIVE  OF  Yd)  SCALFD  BY 

C*  H**J/FACTPRI AL( J) ,  H  THE  CURRFNT  STEP  SIZE. 

C*  ON  THE  FIRST  ENTRY,  THE  CALLER  SUPPLIES 

C*  Ylltl)  ANT  Y(2,I),  UNSCALEO.  (IF  THE  CALL  TO 

C*  DIFSUB  WAS  PRECEDED  BY  A  CALL  TO  PIFMF3,  THE 

C*  USER  NEED  NOT  TOUf"H  Y  AT  ALL).  THE  PROGRAM 

C*  WILL  SCALE  Y(2,I)  BY  H.  ON  ANY  SUBSEQUENT 

C*  ENTRY,  THE  PROGRAM  ASSUMES  THAT  THE  Y  VALUES 

C*  HAVE  NOT  BFEN  CHANGED  SINCE  THE  LAST  EXIT 

C*  FROM  DIFSUB,  AND  WILL  SCALE  THESE  VALUES  IF 

C*  THE  CALLER  HAS  CHANGED  H. 

C*  IF  IT  IS  DESIRED  TO  INTERPOLATE  TO  NON-MESH 

C*  POINTS  THESE  VALUES  CAN  BE  USED.  IF  THE  CURRENT 

C*  STEP  SIZE  IS  H  AND  THE  VALUE  AT  T*E  IS  NEEDED, 

C*  FORM  S  =  E/H  AND  THE  COMPUTE 

C*  NO 

C*  YCIMT+E)  =   SUM   Y(J*1,I  )*S**J 

C*  J*0 

C*  YL       AN  ARRAY  OF  NL  VARIABLES  THAT  APPEAR  ONLY  LINEARLY. 

C*  CALLER  MUST  SUPPLY  VALUES  FOR  THESE  VARIABLES. 

C*  SAVE     AN  ARRAY  OF  LENGTH  AT  LEAST  7*N. 

C*  H        THE  STEP  SIZE  TO  BE  ATTEMPTED  ON  THE  NEXT  STEP. 

C*  H  MAY  BE  AOJUSTED  UP  OR  DOWN  BY  THE  PROGRAM 

C*  IN  ORDER  TO  ACHEIVE  AN  ECONOMICAL  INTEGRATION. 

C*  HOWEVER,  IF  THE  H  PROVIDED  BY  THE  USER  DOES 

C*  NOT  CAUSE  A  LARGER  ERROR  THAN  REQUESTED,  IT 

C*  WILL  BE  USED.   TO  SAVE  COMPUTER  TIME,  THE  USER  IS 

C*  ADVISED  TO  USE  A  FAIRLY  SMALL  STEP  FOR  THE  FIRST 

C*  CALL.   IT  WILL  BE  AUTOMATICALLY  INCREASED  LATFR. 

C*  HMIN     THE  MINIMUM  STEP  SIZE  THAT  WILL  BE  USED  FOR  THE 

C*  INTEGRATION.   NOTE  THAT  ON  STARTING  THIS  MUST  BE 

C*  MUCH  SMALLER  THAN  THE  AVERAGE  H  EXPECTED  SINCE 

C*  A  FIRST  ORDER  METHOD  IS  USED  INITIALLY. 

C*  HMAX     THE  MAXIMUM  ALL3WABLE  STEP  SIZE 

C*  EPS      THE  ERROR  TEST  CONSTANT.   SINGLE  STEP  ERROR  ESTIMATES 

C*  DIVIDED  BY   YMAX(I)   MUST  BE  LESS  THAN  THIS 

C*  IN  THE  EUCLIDEAN  NORM.  THE  STEP  AND/OR  ORDER  IS 

C*  ADJUSTED  TO  ACHEIVE  THIS. 

C*  YMAX     A  VECTOR  OF  LENGTH  NY  WHICH  CONTAINS  THE  MAXIMUM 

C*  OF  EACH  Y  SEEN  SO  FAR.  ON  THE  FIRST  CALL,  THESE 

C*  WILL  BE  INITIALIZED  AS  YMAX(I)  -  MAX( 1 ,  |  Y (  1 ,  I  )  |  ) 

C*  ERROR    A  VECTOR  OF  LENGTH  NY. 

C*  KFLAG    A  COMPLETION  CODE  WITH  THE  FOLLOWING  MEANINGS: 

C*  *l   THE  INTEGRATION  WAS  SUCCESSFUL. 

C*  -1   THE  ERROR  TEST  FAILED  FOR  H  >  HMIN. 

C*  -3   THE  CORRECTOR  FAILED  TO  CONVERGE  FOR 

C*  H  >  HMIN. 

C*  -2       TOO   MANY    FLOATING-POINT    EXCEPTIONS 

C*  OCCURRED    DURING   LAST    STEP. 

C*  -4      THE    CORRECTOR    FAILED    THREE    TIMES    WITH 

C*  EVEN    THE    FIRST-ORDER    METHOD. 


OF  AS 

001 

DFAS 

00? 

DFAS 

003 

DFAS 

004 

OFAS 

005 

DFAS 

OOh 

DFAS 

007 

DFAS 

008 

DFAS 

009 

DFAS 

010 

DFAS 

311 

DFAS 

012 

DFAS 

013 

DFAS 

014 

DFAS 

015 

PCAS 

016 

DFAS 

017 

DFAS 

010 

DFAS 

019 

DFAS 

020 

DFAS 

021 

DFAS 

022 

DFAS 

023 

DFAS 

024 

DFAS 

025 

DFAS 

026 

DFAS 

027 

DFAS 

02  8 

OFAS 

029 

DFAS 

030 

DFAS 

031 

DFAS 

032 

DFAS 

033 

DFAS 

034 

DFAS 

035 

DFAS 

036 

DFAS 

037 

DFAS 

038 

DFAS 

039 

DFAS 

040 

DFAS 

041 

DFAS 

042 

DFAS 

043 

DFAS 

044 

DFAS 

045 

DFAS 

046 

DFAS 

04  7 

DFAS 

048 

DFAS 

049 

DFAS 

050 

DFAS 

051 

DFAS 

052 

DFAS 

053 

DFAS 

054 

DFAS 

055 

DFAS 

056 

DFAS 

057 

DFAS 

05  8 

DFAS 

059 

DFAS 

060 

DFAS 

061 

DFAS 

062 

DFAS 

063 

DFAS 

064 

DFAS 

065 

DFAS 

066 

DFAS 

067 

DFAS 

068 

DFAS 

069 

DFAS 

070 

DFAS 

071 

32 


c* 

r  * 
C* 

r* 
C* 
r  * 

C* 
C* 
C* 
C* 
C* 

c* 

r* 
C* 
C* 
c* 
r* 
c* 
c* 
c* 
c* 


JSTAPT   AN'  INPUT  INCICATOR  WITH  T 

0  PERFORM  THF  FIRS 

INITIALIZES  IT 
AT  IVES  IN  Y<2, 
THF  INTEGRATIO 
ANY  SUBSEQUENT 
WITH  JSTART  = 

1  CONTINUE    FROM    TH 

UNTIL  T  >  TFN'D 
JSTART  IS  SET  TO  NO,  TH 
THE  METHOD,  AT  EXIT. 

MAXPFR  THE  MAXIMUM  DERIVATIVF  TH 
METHOP.  IT  MUST  NOT  EXC 

PW  A  VFCTOR  Oc  LENGTH  N**2*-20  ( 
GENERATED  BY  MATSFT  AND 

YLSV  A     VFCTOR     OF    LENGTH    NL . 

DY  A    VECTOR    "IF    LENGTH    M,     OUT 

ERSV  A    VECTOR    OF    LENGTH    NY. 

Fl  A    VECTOR    OF    LENGTH    N,    OUT 

EQN.VAR  VECTORS    USED    BY    MATSET. 


HE    FOLLOWING    MEANINGS: 
T    STEP.     THF     ROUTINE 
SELF,     SCALES    THF    OER IV- 
I )     AND    THEN    PERFORMS 
N    UNTIL     T    >    TFND. 

CALLS    SHOULD    BE    MADE 
I. 
E    LAST    STEP,     INTEGRATING 

F    CURRENT    ORDER    OF 

AT    SHOULD    BE    USED    IN    THE 
EED    6. 
REAL*4) , 
USFD    BY    MATINV.MATMUL. 

PUT    OF    DIFFUN. 

PUT    OF    MATMUL. 


C**** 

C* 

C**** 


c 

c* 


***************************************************************** 
DIMENSION    T(l),G(ll,Y(7,l),YL(l),SAVE(7,l),YMAX(l) 
DIMENSION    ERROR ( 1 ) , PW  ( 1 )  , YLSV( 1 ) ,DY ( 1 1 , ERSV ( 1 ) 
DIMENSION    Fill  I , EON ( 1),VAP(1),A(7),PEPTST(6,3) 
DATA     PERTST     /4  .  0  ,9.  0  ,  16  .  0,  25  .0  ,  36  .  0  ,  49  .  0  , 9.  0  ,  16  .  0, 

+  25.0,36.0,49.0,64.0,1.0,1.0,0.25, 

♦  2.7889E-2,1.70  569E-3,6.83929E-5/ 
DATA    MF     /?./,    KZILCH    / 0/ 

********************************************* 

THIS  COMMUNICATES  THE  VALUE  OF  N  TO  MATMUL  AND  MATIN3. 

a**********************.**********.********** 

COMMON  N99 

FORMAT  </I5,!4,Dll.2,Dl0.2,6Dl5.5/(31X,6D15.5)) 

FORMAT  (31X.6D15.5) 

FORMAT  ('IMF  =«,I2,«    N  =»iI3i'    NL  =,,I3,'     FPS  =', 

♦  D9.2,'    TFND  =',D9.2,'    H  =',D9.2//I 
FORMAT  (•   NS   NW        H» , 8X , • T( 1 ) • , 8X , 

♦  «Y<  It*)     AND  YL(*)« //) 

NY  =  N-NL 

LENGTH     IN    WORDS    OF     Y,     YL     (FOR    COPY)     : 

LCCPYY    =    NY*14 

LCOPYL    =    NL*2 

N99=N 

IF     (MAXDEP     .GT.    61     MAXDER    =     6 

HINV=2.5 

CALL  DIFFUN  ( T , G, F  I , Y , YL ,HI NV ) 

CALL  MATSET  ( 0, DY , FPS ,EON ,G ,HI NV , 

♦  PW.F1 ,T, VAR,Y,YL) 


0, M,MF,N,NY, 1, 


CALL 

IF  (JST 

WRITE  « 

WRITE  ( 

r  *********** * 

C* 

ON  THE 

c* 

THE  VFC 

r  **  **  *  ....... 

NS  =  0 

NW  =  0 

DO  20  J 

YMAXI 

20 

Y(2,J 

NO  *  1 

BR  =  1. 

ASSIGN 

GO  TO  ? 

r ************ 

C* 

SET  THF 

c* 

TYPF.   1 

c* 

FUP  IS 

c* 

DECRFAS 

r  ************ 

MATIN1      (PW) 
APT    .GT.    01    GO    TO    60 
6,3)     MF.N.NL, EPS, TEND, H 
6,4) 

************************************** 
FIRST    CALL    SET    THE    ORDER    TO    ONE.     INITIALIZE 
TOP     YMAX,    AND     SCALE    THE    FIRST    DERIVATIVES    OF    Y    BY    H. 
************************************** 


=     1,NY 
J)     =    DMAXK l.ODO.DABSIYI  l.JI  )  ) 
)    =     Y(2,J)*H 

0 

100    TO    IRET 

21 

************************************** 

COFPFICIENTS     THAT    DETFRMINF    THE    ORDER    AND    THE    METHOD 
EISA     TEST    FOP    FRRORS    OF    THE    CURRENT    ORDER    NO. 
TO    TEST    FOR     INCREASING    THE    ORDER,     EDWN    FOR 
INC.    THE    ORDER. 
************************************** 


^FAS 

0  72 

HFftS 

073 

DFAS 

074 

OCAS 

075 

DPAS 

0  76 

PFAS 

077 

T=AS 

07R 

HFAS 

079 

PFAS 

OBO 

PPAS 

0R1 

DFAS 

0B2 

P^AS 

0R3 

DFAS 

0R4 

DFAS 

085 

nFAS 

0B6 

DFAS 

087 

DFAS 

088 

DFAS 

089 

DFAS 

090 

DFAS 

091 

DFAS 

092 

**DFAS 

093 

DFAS 

094 

DFAS 

095 

DFAS 

096 

DFAS 

097 

DFAS 

098 

DFAS 

099 

DFAS 

100 

DFAS 

101 

DFAS 

102 

DFAS 

103 

DFAS 

104 

DFAS 

105 

DFAS 

106 

DFAS 

107 

DFAS 

108 

DFAS 

109 

DFAS 

110 

DFAS 

ill 

DFAS 

112 

DFAS 

113 

DFAS 

114 

DFAS 

115 

DFAS 

116 

DFAS 

117 

DFAS 

118 

DFAS 

119 

DFAS 

120 

DFAS 

121 

DFAS 

122 

DFAS 

123 

DFAS 

124 

DFAS 

125 

DFAS 

126 

DFAS 

127 

DFAS 

128 

DFAS 

129 

DFAS 

130 

DFAS 

131 

DFAS 

132 

DFAS 

133 

DFAS 

134 

DFAS 

135 

DFAS 

136 

DFAS 

137 

DFAS 

138 

DFAS 

139 

DFAS 

140 

DFAS 

141 

DFAS 

142 

DFAS 

143 

DFAS 

144 

DFAS 

145 

33 


170   CONTINUE 
221 


226       A(2)    =    -2.45D0 


0<=AS     146 


GOTO    (221,  222,  223,  2  24,225,226),  NO  Dp4S  1*7 

A(2J    =    -l.ODO  DF?-S  1*« 

r,o  to  2?o  °c^  ^9 

222       A(2)     =    -1. 500  j>"S  |5° 

A(3)     =    -C.5D0  °FA^  1S1 

GP    TO    230  nPiS  152 

22*       A(2)     =    -1.83333333333333333  -Ffi<?  l5^ 

A(3)    =    -l.ODO  nPAS  154 

A<4)     =    -0.1666666666666666"'  nFAS  155 

GO    TO    230  nCAS  l5* 

224  A(2)  =  -2.08333333333333233  DFAS  157 
A<3)  =  -1.45833333333333333  DF  AS  15B 
A<4)  =  -0.4166666666666667  DFAS  159 
A(5>  =  -0.04166666666666667  DFAS  160 
GO    Tn    230  nFAS  l61 

225  A(2)  =  -2.8333333333333333  °FAS  l62 
A(3)  =  -1.875D0  DFAS  I63 
A(4)  =  -0.7083333333333333  °FAS  I6* 
A(5)  =  -0.125D0  DFA<>  165 
A(6)  =  -0.0083333333333333333  DFAS  I66 
GO    TO    230  DFAS  167 


DFAS     168 


At3)  =  -2.2555555555555555  nFA5  169 

A(4)  =  -1.2083333333333333  DFAS  I70 

A(5)  =  -0.2430555555555555  OFAS  171 

A(6)  =  -0.02916666666666667  DFAS  172 

A(7)  =  -0.0013888888888888889  DFAS  173 

230   K  =  NO+1  DFAS  17* 

ID0U8  =  NO  nFAS  I75 

FN01  =  0.5/FLDATINQ)  DFAS  l76 

ENQ2  =  0.5/FLDAT(K)  °FAS  I77 

FN03  =  0.5/FL3AT(NQ  ♦  21  DFAS  I78 

PFPSH  =  EPS**2  DFAS  17q 

E    =  PFRTSTINQ, 1) *PEPSH  DFAS  180 

EUP  =  PERTST(NQ,2J*PEPSH  DFAS  181 

FDWN=  PERTST(N0,3I«PEPSH  DFAS  182 

BNO  =  ( EPS*FNQ3 )**2  0FAS  I83 

IWEVAL  =  1  OFAS  184 

GO  TO  TRET,  (100,250,630,751)  OFAS  185 

£***********************•*********•»********•******  DFAS  186 

C*     IF  THE  CALLER  HAS  CHANGED  H,  SCALE  THE  Y  VARIABLES.  DFAS  187 

C*     HNEW  IS  THF  STEP  SIZE  USED  ON  THE  LAST  CALL.  DFAS  188 

£•***********•**********•***********••*************  DFAS  189 

60   IF  (H  .EO.  HNEW)  GO  TO  100  OFAS  190 

R  =  H/HNEW  DFAS  191 

ASSIGN  100  TO  IRET  OFAS  192 

GO  TO  800  DFAS  193 

C**************************************************  DFAS  194 

C*     BEFORE  EXIT,  JSTART  IS  SET  TO  THE  ORDER  OF  THE  METHOD,  DFAS  195 

C*     AND  THE  CURRENT  VALUE  OF  H  IS  SAVFD.  DFAS  196 

£***************•************•***•***************•*  DFAS  197 

70   KFLAG  =  -2  DFAS  198 

80   JSTART  =  NO  DFAS  199 

HNEW  =  H  DFAS  200 

RETURN  DFAS  201 

r I************************************************ ..  DFAS  202 

C*     PRINT  DATA  RELEVANT  TO  THE  STEP,  AND  TAKE  ANOTHER  STEP  DFAS  203 

C*     IF  T  <  TEND.  DFAS  204 

£***«*****•****«************•*»********•********•**  DFAS  205 

90   NS  =  NS*1  DFAS  206 

WRITE  (6,11  NS,NW,H,T (1) ,(Y( 1, I ), 1=1, NY)  DFAS  207 

IF  (NL  .GT.  0)  WRITE  (6,2)  < YL I  I  ) , I  =  1 , NL  )  OFAS  208 

CALL  ANSWER(Y,F1,T)  DFAS  209 

IF    CKNTSPKO)    .GT.    2)    GO    TO    70  DFAS  210 

IF     (KFLAG    .LT.    0)    GO    TO    80  DFAS  211 

IF    (T(l)    .GE.    TEND)    GC   TO    80  DFAS  212 

JSTART    =    1  DFAS  213 

C *********  *****************************************  DFAS  214 

C*     BEGIN  BY  SAVING  INFORMATION  FOR  POSSIBLE  RESTARTS.  DFAS  215 

e*****************************+********************  DFAS  216 

100   CALL  COPYZ  (SAVE,Y,LCOPYY)  DFAS  217 

CALL  COPYZ  (YLSV,YL,LCOPYL)  DFAS  218 
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RACUM  =  1.0  DFAS  219 

KFLAG  =  1  DFAS  220 

HOLD  =  H  DFAS  221 

NOOLD  =  MO  DFAS  22? 

TOLD  =  T(l)  ^FAS  223 

K7ILCH  =  1  DFAS  224 

;*************************************** ********************************DFAS  225 

C*     THIS  SECTION  COMPUTES  THE  PREDICTED  VALUES  BY  EFFECTIVELY  *DPAS  226 

C*     MULTIPLYING  THE  SAVED  INFORMATION  BY  THE  PASCAL  TRIANGLE  *DFAS  227 

C*     MATRIX.  *DFAS  228 

C* *********************************************** ***********************OFAS  229 

250   T(l)  =  TdH-H  DFAS  230 

CALL  S2(T,G)  DFAS  231 

CALL  MATSET  ( 0 , DY , EPS , EON ,G,HI NV ,  0, M,MF , N, N Y ,2  ,  DFAS  232 

+    PW,F1,T, VAR,Y,YL J  DFAS  233 

CALL  MATIN2  <PW)  DFAS  234 

HINV  =  l.DO/H  DFAS  235 

DP  ?60  J  =  2fK  OFAS  236 

J3  =  K*J-1  DFAS  237 

DO  260  Jl  =  JfK  DFAS  238 

J2  =  J3-J1  DFAS  239 

DO  260  I  =  1,NY  DFAS  240 

260          Y(J2,II  =  Y(J2,I)  ♦  YIJ2*l,I)  PpA'S  241 

C******************* ********************************** ******************pFAS  242 

C*       UP  TO  3  CORRECTOR  ITERATIONS  ARE  TAKEN.   CONVFRGENCF  IS  TFSTED  *DFAS  243 

C*       BY  REQUIRING  CHANGES  TO  BE  LESS  THAN  BND  WHICH  IS  DEPENDENT  ON  *DFAS  244 

C*       THE  ERROR  TEST  CONSTANT.  *DFAS  245 

C*           THE  SUM  OF  THE  CCPRECTIONS  IS  ACCUMULATED  IN  THE  ARRAY  *DFAS  246 

C*       ERROR(I).   IT  IS  EQUAL  TO  THE  K-TH  DERIVATIVE  OF  Y  MULTIPLIED  *DFAS  247 

C*       BY   H**K/(FACTORI AL  (K-1)*A(K) ),  AND  IS  THERE  =  ORE  PROPORTIONAL  *DFAS  248 

C*       TO  THE  ACTUAL  ERRORS  TO  THE  LOWEST  POWER  OF  H  PRESENT.  (H«*K)  *DFAS  249 
£**** ******** ******* ****************************************************  op  as  250 

DO  270  I  =  1,NY  DFAS  251 

270     ERROR(I)  =  0.0  DFAS  252 

DO  430  L  =  1,3  DFAS  253 

CALL  DIFFUN  d , G, DY  ,Y , YL, HI NV I  DFAS  254 

C ***** ******** **********************************************************  (if AS  255 

C*       IF  THERE  HAS  BEEN  A  CHANGE  OF  ORDER  OR  THERE  HAS  BEFN  TROUBLE  *OFAS  256 

C*       WITH  CONVERGENCE,  PW  IS  RE-EVALUATED  PRIOR  TO  STARTING  THE  DFAS  257 

C*       CORRECTOR  ITERATION  IN  THE  CASE  OF  STIFF  METHODS.   IWEVAL  IS  *DFAS  258 

C*       THEN  SET  TO  -1  AS  AN  INDICATOR  THAT  IT  HAS  BEEN  DONE.  *DFAS  259 

C ********* ******** *********** *******  ************************** **********qfaS  260 

IF  (IWEVAL  .LT.  1  .AND.  L  .GT.  I)  GO  TO  280  DFAS  261 

CALL  MATSET  ( A ( 2) , CY , EPS , FQN, G,H I NV ,  0, M , MF , N.NY , 3 ,  DFAS  262 

♦    PW,F1,T,VAR,Y,YL>  DFAS  263 

CALL  MATIN3  IPW)  DFAS  264 

KFLAG  *  1  DFAS  265 

IWEVAL  «  -1  DFAS  266 

NW  ■  NW*1  DFAS  267 

280       CALL  MATMUL  «PW,DY,F1)  OFAS  268 

IF  (NL  .LE.  0)  GO  TO  300  DFAS  269 

DO  290  I  -  l.NL  DFAS  270 

290       YL(I)  -  YL(II  -  FlU+NY)  DFAS  271 

300     CONTINUE  DFAS  272 

DEL  -  0.000  DFAS  273 

DO  420  I  -  I, NY  DFAS  274 

Yd, I)  »  Yd, I)  -  F1III  OFAS  275 

Y(2,I»  -  Y(2,I»  ♦  A<2»*F1(I)  DFAS  276 

ERROR(I)  -  ERROR(I)  ♦  F1(I»  DFAS  277 

DEL  -  DEL  ♦  (Fill )/YMAX< I) )**2  DFAS  278 

420       CONTINUE  DFAS  279 

IFIL.GE.2)  BR  -  DM  AX  1 1  .9*BR , DEL/DEL  1  I  DFAS  280 

DELI  ■  DEL  DFAS  281 

IF(DMIN1IDEL,BR*0EL*2.0).LE.BND|  GO  TO  490  DFAS  282 

430     CONTINUE  OFAS  283 

C************* ******************************** **************************DFAS  284 

C*            THE    CORRECTOR    ITERATION    FAILEO    TO    CONVERGE     IN    3    TRIES.       VARIOUS  *DFAS    285 

C*            POSSIBILITIES    ARE    CHFCKED    FOR.       IF    H    IS     ALREADY    HMIN    AND  *DFAS    286 

C*            THIS     IS    EITHER    ADAMS    METHOD    OR    THE    STIFF    METHOD    IN    WHICH    THE  *0FAS    287 

C*            MATRIX    PW    HAS    ALREADY    BEEN    RE-EVALUATED,     A    NO    CONVERGENCE    EXIT  *DFAS    288 

C*            IS    TAKEN.    OTHERWISE    THE    MATRIX    PW    IS    RE-EVALUATED    AND/OR    THE  *DFAS    289 

C*            STEP    IS    RFDUCED    TO    TRY    AND    GET    CONVERGENCE.  *DFAS    290 

C************* **************************** ******************************[)FAS     291 
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440       T(l)    =    TOLD  nc4S  292 

IF     (IWEVAL)     445,455,450  OFAS  293 

445        IF     (H     .LF.     HMI N* 1 . 000000 1 )     GO    TO    460  OFAS  2°4 

450       RACUM    =     PACUM*0.25  OFAS  295 

455       CONTINUE  ^FAS  29f 

GO    TO    750  OFAS  277 

460       KFLAG    =     -3  OFtS  29R 

470      CALL    COPYZ     ( Y, SAVF , LC CPYY )  ~>FAS  2C9 

CALL    COPYZ     (YL.YLSV.LCOPYL)  OFAS  300 

H    =    HOLD  DFAS  301 

NO  =  NQOLD  DFAS  302 

GO  TO  90  nF:AS  303 

r *************************************** ***********  OFAS  304 

C*     THF  CORRECTOR  CONVFRGEC  ANO  NOW  THF  ERROR  TEST  IS  MADE.             OFAS  305 

C **************************************************  OFAS  306 

490   P  =  0.0  OFAS  307 

00  500  I  =  1,M1  OFAS  308 

500     0=0*  (ERROR (I )/YMAX< I ) )**2  OFAS  309 

IWEVAL  =  0  OFAS  310 

IF  (D  .GT.  E)  GO  TO  540  OFAS  311 

C **************************************************  DCAS  312 

C*     THE  ERROR  TEST  IS  OKAY,  SO  THE  STEP  IS  ACCEPTED.  IF  IDOUR  OFAS  313 

C*     NOW  BECOMES  NEGATIVE,  A  TEST  IS  MADE  TO  SEE  IF  THE  STFP             OFAS  314 

C*     SIZE  CAN  BE  INCREASED  AT  THIS  ORDFP  OP  ONE  HIGHER  OR  OFAS  315 

C*     LOWER.  THE  CHANGE  IS  MADE  ONLY  IF  THE  STFP  CAN  BE  IN-               OFAS  316 

C*     CREASED  RY  AT  LEAST  10*.  IDOUB  IS  SET  TO  NO  TO  PREVENT              DFAS  317 

C*     FURTHER  TFSTING  FOR  A  WHILE.  IF  NO  CHANGE  IS  MADE,  IDOUR  DFAS  318 

C*      IS  SFT  TO  9.  OFAS  319 

C*************************************+************  DFAS  320 

IF     (K    .LT.     3)     GO    TO    5  20  DFAS  321 

DO    510    J    =    3,K  DFAS  322 

DO    510    I     =    1,NY  DFAS  323 

510                   Y(J,I)     =    Y(J,I)     ♦    A(J)*ERRnR(I)  DFAS  324 

520      KFLAG    =    >1  DFAS  325 

IDOUB    =     IDOUB-1  DFAS  326 

IF    (IDOUB)     550,525,700  DFAS  327 

525       DO    530    I    =     1.M1  OFAS  328 

530            ERSVdl    =    ERROR(I)  DFAS  329 

GO    TO    700  OFAS  330 

C**************************************************  DFAS  331 

C*     THE  ERROR  TEST  FAILED.  IF  JSTART  =  0,  THE  DERIVATIVES  DFAS  332 

C*     IN  THE  SAVE  ARRAY  ARE  UPDATED.  TESTS  ARE  THEN  MADE  TO  DFAS  333 

C*     FIX  THE  STEP  SIZE  AND  PERHAPS  REDUCE  THE  ORDER.  AFTER  OFAS  334 

C*     RESTORING  AND  SCALING  THE  Y  VARIABLES,  THE  STEP  IS  OFAS  335 

C*     RFTRIED.  DFAS  336 

C**************************************************  DFAS  337 

540   IF  (JSTART  .GT.  OJ  GO  TO  548  DFAS  338 

DO  544  I  =  l.NY  DFAS  339 

544     SAVE(2,l)  =  Y(2,I)  DFAS  340 

548   KFLAG  «  KFLAG  -  2  DFAS  341 

IF  (H  .LE.  HMINI  GO  TO  740  DFAS  342 

T(l)  »  TOLD  OFAS  343 

IF  (KFLAG  .LE.  -5)  GO  TO  720  OFAS  344 

550   PR2  *  (D/E)**EN02*1.2  DFAS  345 

L  «  0  DFAS  346 

IF  (NO  .LE.  1)  GO  TO  570  OFAS  347 

D  -  0  DFAS  348 

DO  560  J  ■  1,M1  DFAS  349 

560       D  ■  D*( Y(K,J)/YMAX(J>)**2  DFAS  350 

PR1  -  (D/E0WN)«*EN01*1.3  DFAS  351 

IF  (PR1  .GE.  PR2)  GC  TO  570  DFAS  352 

PR2  *  PR1  DFAS  353 

L  -  -1  OFAS  354 

570   IF  (KFLAG  .LT.  0   .OR.   NO  . GE .  MAXDERI  GO  TO  590  DFAS  355 

D  =  0  OFAS  356 

DO  580  J  *  1,M1  DFAS  357 

580       D  -  DM (ERROR(J)-ERSV( J) )/YMAX(J) )**2  DFAS  358 

PR1  =  (D/EUPI**EN03*1.4  DFAS  359 

IF  (PRl  .GE.  PR2)  GC  TO  590  DFAS  360 

PR2  -  PRl  DFAS  361 

L  *  1  DFAS  362 

590   R  =  1.0/AMAXHPR2,  l.E-5)  DFAS  363 

IF  (KFLAG  .LT.  0   .OR.   R  .GE.  1.1)  GO  TO  600  DFAS  364 

IDOUB  -  9  DFAS  365 

GO  TO  700  OFAS  366 


36 


600       NEWO    =    NO+L  DFAS    367 

K    =    NFWQ+1  DFAS     368 

IF     (NFWO    .LE.    NQI    GO    TO    620  DFAS     369 

Rl    =    A(NFWQ)/FLDAT(NEV<Q)  DFAS    370 

DO    610    J    =     1,NY  DFAS     371 

610  Y(K,J>    =    ERR0R(J)*P1  DFAS     37? 

620       C"NTINUF  DFAS     373 

C***#*********#*****************************+******  DFAS    374 

C*     IF  THE  STEP  WAS  OKAY,  SCALE  THE  Y  VARIABLES  IN  ACCORDANCE  DFAS  375 

r*     WITH  THE  NEW  VALUE  OF  H.  IF  KFLAG  <  0,  HOWEVER,  USE  THE  DFAS  376 

C*     SAVED  VALUES  (IN  SAVE  AND  YLSV).  IN  EITHER  CASE,  IF  THE  DFAS  377 

C*     ORDER  HAS  CHANGED  IT  IS  NECESSARY  TO  FIX  CERTAIN  PARAMETFRS  DFAS  378 

C*     BY  CALLING  THE  PROGRAM  SEGMENT  AT  LINE  170.  DFAS  379 

r  ****  «•****»#**  *****.*,*****************.******<.****:.!<**.  DFAS  380 

IDOUB  =  NO  DFAS  381 

IF  (NEWO  .EQ.  NO)  GO  TO  630  DFAS  382 

NO  =  NFWO  DCAS  383 

ASSIGN  630  TO  IPET  DFAS  384 

GO  TO  170  DFAS  385 

630   IF  (KFLAG  .GT.  0)  GO  TO  670  DFAS  386 

RACUM  =  RACUM*R  DFAS  387 
GO  TO  750                                                          '     DFAS  388 

670   P  =  DMAX1(DMIN1(HMAX/H,R),HMIN/H)  DFAS  389 

H  =  H*R  DFAS  390 

IWEVAL  =  1  DFAS  391 

ASSIGN  700  TO  IRET  DFAS  392 

GO  TO  800  DFAS  393 

700   DO  710  I  =  1,M1  DFAS  394 

710     YMAX(I)  =  DMAXK  YMAXU  )  ,DARS(Y(  1,  I  )  ))  DFAS  395 

GO  TO  90  DFAS  396 

C  *********************************************  *****  DFAS  397 

C*     THE  ERROR  TEST  HAS  NOW  FAILED  THREE  TIMES,  SO  THE  DFAS  398 

C*     DERIVATIVES  ARE  IN  BAD  SHAPE.  RETUPN  TO  FIRST-ORDER  DFAS  399 

C*     METHOD  AND  TRY  AGAIN.  OF  COURSE,  IF  NQ  =  1  ALREAOY,  DFAS  400 

C*     THEN  THERE  IS  NO  HOPE  AND  WE  EXIT  WITH  KFLAG  =  -4.  DFAS  401 

C**************************************************  DFAS  402 

720   IF  (NO  .EO.  I)  GO  TO  735  DFAS  403 

NO  =  1  DFAS  404 

IDOUB  =  1  DFAS  405 

ASSIGN  751  TO  IPET  DFAS  406 

GO  TO  170  DFAS  407 

735   NOOLD  *  1  DFAS  408 

KFLAG  «  -4  DFAS  409 

GO  TO  470  DFAS  410 

740   KFLAG  =  -1  DFAS  411 

GO  TO  90  DFAS  412 

r *»********•**•******»**•••*•*•**•.»•**»..».* .... »*  DFAS  413 

C*     THIS  SECTION  RESTORES  THE  SAVED  VALUES  OF  Y  AND  YL ,  SCALING  DFAS  414 

C*     THE  Y  DERIVATIVES  AS  NECESSARY,  AND  THEN  RETURNS  TO  THE  DFAS  415 

C*     PREDICTER  LOOP.  DFAS  416 

r ****+**********♦*♦******* *•*•*•«*•******* *********  DFAS  417 

750  H  ■  HOLD*RACUM  DFAS  418 
H    ■    DMAX1(HMIN,0MIN1(H,HMAX) )  DFAS    419 

751  RACUM  -  H/HOLO  DFAS  420 
Rl  «  1.0  DFAS  421 
DO    760    J    ■    2,K  DFAS    422 

Rl    «    R1*RACUM  DFAS    423 

00    760    I    ■    1,NY  DFAS    424 

760  Y(J,IJ    ■    SAVE(J,II*Rl  DFAS    425 

DO    770    I    -    liNV  DFAS    426 

770  Yd, I)    -    SAVE(l.I)  DFAS    427 

CALL    COPYZ     (YL.YLSV.LCOPYll  DFAS    428 

IWEVAL    -    1  DFAS    429 

GO    TO    250  DFAS    430 

r I*.**.*.****.***.***.**..*.*.*.,,*.*.*****. ****...  DFAS    431 

C*     THIS  SECTION  SCALES  THE  Y  DERIVATIVES  BY  R.  DFAS  432 

r *********************************** ***************  DFAS  433 

800   Rl  «  1.0  DFAS  434 

DO  810  J  -  2,K  DFAS  435 

PI  *  R1*R  DFAS  436 

DO  810  I  »  1,NY  OFAS  437 

810       Y(J,I)  «  Y(J,I)*R1  DFAS  438 

GO  TO  IRET,  (100,700»  DFAS  439 
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C* 

r ****** 

ENTRY    RFDSUB 


************ *****, 


♦•♦•♦A************************************** 

PNTFR    HERE    TO    RESET    VALUES     IN    PREP    FOP    AN    AUTO-RESTAPT . 
I******************************************* 

ENTRY    RFDSUB 

IF     INTERRUPT    OCCURPFD    IN    MATIN1,     Nn    RFSTORATIDN    NECESSARY: 

IF    (KZILCH    .EO.    0)    RETURN 


"0    910    J    =     1,NY 

Y(  I, J)    =    SAVEIl, J> 
910  Y(2,J)     =    SAVEI2, JJ/HOLD 

CALL    CPPYZ     (YL,YLSV,LCOPYL) 
T(  1)     =    TOLD 
RETURN 
'      END 
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